Exploración de los datos
Librerías
# Load necessary libraries
library(arm) # Regression modeling
## Loading required package: MASS
## Loading required package: Matrix
## Loading required package: lme4
##
## arm (Version 1.14-4, built: 2024-4-1)
## Working directory is /Users/celeste.solorzano/Desktop/Maestría Estadistica/2024/II Ciclo/Modelo Mixtos/Proyecto final
library(corrplot) # Correlation plots
## corrplot 0.95 loaded
##
## Attaching package: 'corrplot'
## The following object is masked from 'package:arm':
##
## corrplot
library(dplyr) # Data manipulation
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2) # Data visualization
library(readxl) # Read Excel files
library(tidyverse) # Collection of data science tools
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ✔ readr 2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ dplyr::select() masks MASS::select()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidyr) # Tidying data
Datos
# Read data
data0 <- read_excel("~/Desktop/Maestría Estadistica/2024/II Ciclo/Modelo Mixtos/Proyecto final/data.xlsx")
# Select relevant data columns
data = data0[,-c(1,3:8)]
# Convert variables to factors
data$Parte = as.factor(data$Parte)
data$Epoca = as.factor(data$Epoca)
data$Clase = as.factor(data$Clase)
# Aggregate 'Parte' into 'Seccion'
data$Seccion <- factor(
ifelse(data$Parte %in% c("I1", "I2", "I3"), "Alta",
ifelse(data$Parte %in% c("M1", "M2", "M3"), "Media", "Baja"))
)
# Aggregate 'Clase' into 'Conta'
data$Conta <- (
ifelse(data$Clase %in% c("1", "2"), 0,
ifelse(data$Clase %in% c("3", "4", "5"), 1, 0))
)
Exploración
# Descriptive analysis of the dependent variable
tabla_frecuencias <- table(data$ICA)
tabla_proporciones <- prop.table(tabla_frecuencias)
barplot(table(data$ICA), main="", xlab="Valor", ylab="Frecuencia", col="#A9C6E7")
# Correlation plot
corrplot(cor(data[,-c(1,2,28,29,30)]), method = "square", type = "lower")
# Convert ICA to factor and get the number of unique categories
data$ICA = as.factor(data$ICA)
num_categories <- length(unique(data$ICA))
# Create a color palette for the categories
formal_blue_palette <- colorRampPalette(c("#A9C6E7", "#6FA3F0", "#004B87", "#002A5C"))(num_categories)
# Function to create a bar plot using tidy evaluation
create_bar_plot <- function(x_var, x_label) {
ggplot(data, aes(x = .data[[x_var]], fill = ICA)) + # Use .data for tidy evaluation
geom_bar(position = "fill", color = "black") +
labs(title = "", x = x_label, y = "Proporción") +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = formal_blue_palette) +
theme_minimal() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
}
# Generate plots for each variable
plot_parte <- create_bar_plot("Parte", "Sección de la Cuenca")
plot_epoca <- create_bar_plot("Epoca", "Época")
plot_seccion <- create_bar_plot("Seccion", "Sección de la Cuenca")
# Display plots
print(plot_parte)
print(plot_epoca)
print(plot_seccion)
# Convert response variable to numeric
data$ICA = as.numeric(data$ICA)
# Standardize numeric variables
datos_est <- data %>% mutate(across(where(is.numeric) & (3:27), scale))
# Select a sample of variables
datos_grap = datos_est[, c(1:2, 8, 22, 27:30)]
# Convert the dataframe to long format for ggplot
datos_long <- datos_grap %>%
pivot_longer(cols = 3:4) %>%
janitor::clean_names() # Clean column names
# Function to create correlation plots
create_correlation_plot <- function(color_var, title) {
ggplot(datos_long, aes(x = value, y = ica, color = !!sym(color_var))) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE) +
facet_wrap(~ name) +
labs(title = title) +
theme_minimal()
}
# Generate correlation plots for different color variables
plot_climate <- create_correlation_plot("clase", "Correlación de Variables con Calidad de Agua por Clase")
plot_interpretation <- create_correlation_plot("interpretacion_de_calidad", "Correlación de Variables con Calidad de Agua por Interpretación de Calidad")
plot_epoca <- create_correlation_plot("epoca", "Correlación de Variables con Calidad de Agua por Época Climática")
plot_parte <- create_correlation_plot("parte", "Correlación de Variables con Calidad de Agua por Parte de la Cuenca")
plot_seccion <- create_correlation_plot("seccion", "Correlación de Variables con Calidad de Agua por Sección de la Cuenca")
# Display plots
print(plot_climate)
## `geom_smooth()` using formula = 'y ~ x'
print(plot_interpretation)
## `geom_smooth()` using formula = 'y ~ x'
print(plot_epoca)
## `geom_smooth()` using formula = 'y ~ x'
print(plot_parte)
## `geom_smooth()` using formula = 'y ~ x'
print(plot_seccion)
## `geom_smooth()` using formula = 'y ~ x'
# Function to create heatmap
create_heatmap <- function(x_var, title) {
ggplot(datos_long, aes_string(x = x_var, y = "epoca", fill = "ica")) +
geom_tile() +
scale_fill_gradient(low = "lightblue", high = "darkblue") +
labs(title = title,
x = "Sección de la Cuenca",
y = "Época",
fill = "ICA") +
theme_minimal()
}
# Generate heatmaps for 'parte' and 'seccion'
heatmap_parte <- create_heatmap("parte", "Heatmap del Índice de Calidad de Agua por Parte")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
heatmap_seccion <- create_heatmap("seccion", "Heatmap del Índice de Calidad de Agua por Sección")
# Display heatmaps
print(heatmap_parte)
print(heatmap_seccion)
Modelos
Librerías
# Load necessary libraries
library(MASS)
library(MCMCglmm)
## Loading required package: coda
##
## Attaching package: 'coda'
## The following object is masked from 'package:arm':
##
## traceplot
## Loading required package: ape
##
## Attaching package: 'ape'
## The following object is masked from 'package:dplyr':
##
## where
## The following object is masked from 'package:arm':
##
## balance
library(R2jags)
## Loading required package: rjags
## Linked to JAGS 4.3.2
## Loaded modules: basemod,bugs
##
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
##
## traceplot
## The following object is masked from 'package:arm':
##
## traceplot
library(dplyr)
library(leaps)
library(readxl)
library(sp)
load.module("dic")
## module dic loaded
load.module("glm")
## module glm loaded
Selección de variables
# Scale selected columns (3 to 26) in the dataset
data_aux <- data %>% mutate(across(c(3:26), ~ scale(., center = TRUE, scale = FALSE)))
# Remove specific columns from the scaled dataset
data_aux1 <- data_aux[, -c(28, 29)]
# Exclude the "ICA" variable from the dataset
variables <- data_aux1[, -which(names(data_aux1) == "ICA")]
# Define the outcome variable as "ICA"
outcome <- data_aux1$ICA
# Perform best subset selection for regression
out_1 <- regsubsets(ICA ~ ., data = data_aux1)
## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 3 linear dependencies found
## Reordering variables and trying again:
# Output the summary of the best subset selection
summary(out_1)
## Subset selection object
## Call: regsubsets.formula(ICA ~ ., data = data_aux1)
## 37 Variables (and intercept)
## Forced in Forced out
## ParteF2 FALSE FALSE
## ParteF3 FALSE FALSE
## ParteI1 FALSE FALSE
## ParteI2 FALSE FALSE
## ParteI3 FALSE FALSE
## ParteM1 FALSE FALSE
## ParteM2 FALSE FALSE
## ParteM3 FALSE FALSE
## EpocaSeca FALSE FALSE
## EpocaTransicion FALSE FALSE
## `Caudal m3/s` FALSE FALSE
## `T aire °C` FALSE FALSE
## `Vel viento m/s` FALSE FALSE
## `Hum Rel %` FALSE FALSE
## pH FALSE FALSE
## `Temp agua (°C)` FALSE FALSE
## `OD (mg/l)` FALSE FALSE
## `PSO (%)` FALSE FALSE
## `DBO Sol (mgO2/l)` FALSE FALSE
## `DBO Total (mgO2/l)` FALSE FALSE
## `DBO rapida` FALSE FALSE
## `NO3 (ug/l)` FALSE FALSE
## `NO2 (ug/l)` FALSE FALSE
## `SST (mg/l)` FALSE FALSE
## `Cond (µS/cm)` FALSE FALSE
## `NH4 (ugNH4/l)` FALSE FALSE
## `Alca (mg/l)` FALSE FALSE
## `P Tot (mgPO/l)` FALSE FALSE
## `P sol (mgPO/l)` FALSE FALSE
## `N Tot (ugN/l)` FALSE FALSE
## `N Org (ugN/l)` FALSE FALSE
## `Prom precip (mm)` FALSE FALSE
## `Precip acum men (mm)` FALSE FALSE
## Conta FALSE FALSE
## `DBO Lenta` FALSE FALSE
## SeccionBaja FALSE FALSE
## SeccionMedia FALSE FALSE
## 1 subsets of each size up to 9
## Selection Algorithm: exhaustive
## ParteF2 ParteF3 ParteI1 ParteI2 ParteI3 ParteM1 ParteM2 ParteM3
## 1 ( 1 ) " " " " " " " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " " " " " " "*"
## 5 ( 1 ) " " " " " " " " " " " " " " " "
## 6 ( 1 ) " " " " " " " " " " " " " " "*"
## 7 ( 1 ) " " " " " " " " " " " " " " "*"
## 8 ( 1 ) " " " " " " " " " " " " " " "*"
## 9 ( 1 ) " " " " " " " " " " " " "*" "*"
## EpocaSeca EpocaTransicion `Caudal m3/s` `T aire °C` `Vel viento m/s`
## 1 ( 1 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " "
## 5 ( 1 ) " " " " " " " " " "
## 6 ( 1 ) " " " " " " " " " "
## 7 ( 1 ) " " " " " " " " "*"
## 8 ( 1 ) " " " " " " " " "*"
## 9 ( 1 ) " " " " " " "*" "*"
## `Hum Rel %` pH `Temp agua (°C)` `OD (mg/l)` `PSO (%)`
## 1 ( 1 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " "*"
## 3 ( 1 ) " " " " " " " " "*"
## 4 ( 1 ) " " " " " " " " "*"
## 5 ( 1 ) "*" " " " " " " "*"
## 6 ( 1 ) "*" " " " " " " "*"
## 7 ( 1 ) "*" " " " " " " "*"
## 8 ( 1 ) "*" " " "*" " " "*"
## 9 ( 1 ) "*" " " " " " " "*"
## `DBO Sol (mgO2/l)` `DBO Total (mgO2/l)` `DBO rapida` `DBO Lenta`
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) "*" " " " " " "
## 4 ( 1 ) "*" " " " " " "
## 5 ( 1 ) "*" " " " " " "
## 6 ( 1 ) "*" " " " " " "
## 7 ( 1 ) "*" " " " " " "
## 8 ( 1 ) "*" " " " " " "
## 9 ( 1 ) "*" " " " " " "
## `NO3 (ug/l)` `NO2 (ug/l)` `SST (mg/l)` `Cond (µS/cm)` `NH4 (ugNH4/l)`
## 1 ( 1 ) " " " " " " " " " "
## 2 ( 1 ) " " " " " " " " " "
## 3 ( 1 ) " " " " " " " " " "
## 4 ( 1 ) " " " " " " " " " "
## 5 ( 1 ) " " " " " " " " "*"
## 6 ( 1 ) " " " " " " " " "*"
## 7 ( 1 ) " " " " " " " " "*"
## 8 ( 1 ) " " " " " " " " "*"
## 9 ( 1 ) " " " " " " " " "*"
## `Alca (mg/l)` `P Tot (mgPO/l)` `P sol (mgPO/l)` `N Tot (ugN/l)`
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " " " " " " "
## 4 ( 1 ) " " " " " " " "
## 5 ( 1 ) " " " " " " " "
## 6 ( 1 ) " " " " " " " "
## 7 ( 1 ) " " " " " " " "
## 8 ( 1 ) " " " " " " " "
## 9 ( 1 ) " " " " " " " "
## `N Org (ugN/l)` `Prom precip (mm)` `Precip acum men (mm)` SeccionBaja
## 1 ( 1 ) " " " " " " " "
## 2 ( 1 ) " " " " " " " "
## 3 ( 1 ) " " " " " " " "
## 4 ( 1 ) " " " " " " " "
## 5 ( 1 ) " " " " " " " "
## 6 ( 1 ) " " " " " " " "
## 7 ( 1 ) " " " " " " " "
## 8 ( 1 ) " " " " " " " "
## 9 ( 1 ) " " " " " " " "
## SeccionMedia Conta
## 1 ( 1 ) " " "*"
## 2 ( 1 ) " " "*"
## 3 ( 1 ) " " "*"
## 4 ( 1 ) " " "*"
## 5 ( 1 ) " " "*"
## 6 ( 1 ) " " "*"
## 7 ( 1 ) " " "*"
## 8 ( 1 ) " " "*"
## 9 ( 1 ) " " "*"
# Plot the results of the best subset selection
plot(out_1)
# Select variables using stepwise AIC method on a pooled linear model
modelo_seleccionado <- stepAIC(lm(ICA ~., data = data_aux1), direction = "both", trace = FALSE)
# Output the summary of the selected model
summary(modelo_seleccionado)
##
## Call:
## lm(formula = ICA ~ Parte + `Caudal m3/s` + `Vel viento m/s` +
## `Hum Rel %` + pH + `OD (mg/l)` + `PSO (%)` + `DBO Total (mgO2/l)` +
## `DBO rapida` + `NO3 (ug/l)` + `Cond (µS/cm)` + `NH4 (ugNH4/l)` +
## `P Tot (mgPO/l)` + `N Tot (ugN/l)` + `N Org (ugN/l)` + `Precip acum men (mm)` +
## Conta, data = data_aux1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75441 -0.34290 0.01357 0.33878 0.82511
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.049e+00 4.475e-01 9.048 6.07e-10 ***
## ParteF2 5.806e-01 3.529e-01 1.645 0.110702
## ParteF3 4.416e-02 3.873e-01 0.114 0.910019
## ParteI1 1.586e-01 6.138e-01 0.258 0.797878
## ParteI2 -1.247e-01 5.836e-01 -0.214 0.832297
## ParteI3 3.242e-01 5.184e-01 0.625 0.536571
## ParteM1 2.658e-01 4.157e-01 0.639 0.527558
## ParteM2 -7.421e-02 4.385e-01 -0.169 0.866804
## ParteM3 -6.608e-01 5.122e-01 -1.290 0.207210
## `Caudal m3/s` -4.064e-01 1.597e-01 -2.545 0.016489 *
## `Vel viento m/s` 3.041e-01 9.254e-02 3.286 0.002663 **
## `Hum Rel %` 4.707e-02 1.665e-02 2.828 0.008410 **
## pH 3.823e-01 2.364e-01 1.617 0.116769
## `OD (mg/l)` -1.774e-01 1.416e-01 -1.253 0.220226
## `PSO (%)` -2.081e-02 1.046e-02 -1.990 0.056044 .
## `DBO Total (mgO2/l)` -4.471e-02 1.137e-02 -3.934 0.000479 ***
## `DBO rapida` 7.819e-02 1.380e-02 5.665 4.01e-06 ***
## `NO3 (ug/l)` 1.174e-04 8.781e-05 1.337 0.191747
## `Cond (µS/cm)` -6.942e-03 2.566e-03 -2.705 0.011312 *
## `NH4 (ugNH4/l)` 2.076e-04 4.422e-05 4.695 5.91e-05 ***
## `P Tot (mgPO/l)` 2.424e-01 1.590e-01 1.525 0.138158
## `N Tot (ugN/l)` -4.267e-04 2.976e-04 -1.434 0.162230
## `N Org (ugN/l)` 4.092e-04 3.010e-04 1.360 0.184430
## `Precip acum men (mm)` 1.071e-03 1.018e-03 1.051 0.301759
## Conta 2.422e+00 3.348e-01 7.234 5.77e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5459 on 29 degrees of freedom
## Multiple R-squared: 0.9694, Adjusted R-squared: 0.9442
## F-statistic: 38.34 on 24 and 29 DF, p-value: 4.624e-16
Modelo 1: Pooling
# Define and configure the JAGS model
jags_model_code <- "
model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i], tau.y) # Likelihood: y follows a normal distribution
mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i] # Linear predictor
e.y[i] <- y[i] - mu[i] # Residuals
}
beta0 ~ dnorm(0, 0.001) # Prior for intercept
beta1 ~ dnorm(0, 0.001) # Prior for slope of x1
beta2 ~ dnorm(0, 0.001) # Prior for slope of x2
tau.y <- pow(sigma.y, -2) # Precision parameter
sigma.y ~ dunif(0, 100) # Prior for standard deviation
}"
writeLines(jags_model_code, con = "modelo_jags.bug") # Save model code to file
# Prepare data for JAGS
data_jags <- list(
y = data$ICA, # Response variable
x1 = as.numeric(data$`Temp agua (°C)`), # Predictor variable 1
x2 = as.numeric(data$`P sol (mgPO/l)`), # Predictor variable 2
n = nrow(data) # Number of observations
)
# Run the JAGS model
jags_model_results <- jags(
data = data_jags,
inits = function() list(beta0 = rnorm(1), beta1 = rnorm(1), beta2 = rnorm(1), sigma.y = runif(1, 0, 100)), # Initial values
parameters.to.save = c("beta0", "beta1", "beta2", "sigma.y", "e.y"), # Parameters to monitor
model.file = "modelo_jags.bug", # Model file
n.iter = 3000, # Total iterations
n.burnin = 1000, # Burn-in period
n.thin = 2, # Thinning interval
n.chains = 3 # Number of chains
)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 54
## Unobserved stochastic nodes: 4
## Total graph size: 344
##
## Initializing model
print(jags_model_results) # Print results
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
## 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 2
## n.sims = 3000 iterations saved. Running time = 0.176 secs
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta0 -8.444 2.960 -14.278 -10.368 -8.480 -6.525 -2.743 1.002 1700
## beta1 0.600 0.137 0.332 0.512 0.604 0.689 0.874 1.001 2700
## beta2 1.372 0.593 0.186 0.986 1.369 1.760 2.552 1.003 820
## e.y[1] -1.687 0.696 -3.047 -2.141 -1.675 -1.232 -0.324 1.002 1600
## e.y[2] -1.082 0.622 -2.285 -1.486 -1.071 -0.671 0.122 1.002 1600
## e.y[3] -1.516 0.667 -2.830 -1.951 -1.519 -1.078 -0.210 1.002 1700
## e.y[4] -2.416 0.486 -3.356 -2.743 -2.419 -2.094 -1.467 1.002 1700
## e.y[5] -3.836 0.410 -4.632 -4.109 -3.841 -3.565 -3.034 1.002 1800
## e.y[6] -5.336 0.270 -5.858 -5.513 -5.344 -5.155 -4.796 1.001 3000
## e.y[7] -1.535 0.407 -2.357 -1.801 -1.531 -1.259 -0.737 1.001 3000
## e.y[8] 1.065 0.317 0.429 0.853 1.062 1.282 1.682 1.001 3000
## e.y[9] 1.765 0.358 1.047 1.532 1.766 2.008 2.466 1.001 3000
## e.y[10] 3.864 0.361 3.170 3.626 3.858 4.101 4.572 1.001 2100
## e.y[11] 4.044 0.389 3.296 3.783 4.036 4.302 4.805 1.001 2100
## e.y[12] 1.864 0.361 1.170 1.626 1.858 2.101 2.572 1.001 2500
## e.y[13] 2.084 0.276 1.550 1.898 2.080 2.269 2.638 1.001 3000
## e.y[14] -0.036 0.271 -0.565 -0.217 -0.043 0.147 0.504 1.001 3000
## e.y[15] -1.156 0.268 -1.678 -1.336 -1.164 -0.973 -0.613 1.001 3000
## e.y[16] 0.045 0.483 -0.919 -0.272 0.047 0.371 0.995 1.001 3000
## e.y[17] -1.255 0.541 -2.366 -1.610 -1.258 -0.890 -0.197 1.001 3000
## e.y[18] -0.435 0.577 -1.620 -0.817 -0.440 -0.044 0.690 1.001 3000
## e.y[19] 0.924 0.370 0.214 0.679 0.917 1.166 1.647 1.002 1900
## e.y[20] 0.624 0.327 -0.015 0.410 0.615 0.843 1.267 1.001 2100
## e.y[21] 0.864 0.361 0.170 0.626 0.858 1.101 1.572 1.002 1900
## e.y[22] -0.136 0.361 -0.830 -0.374 -0.142 0.101 0.572 1.002 1900
## e.y[23] -0.496 0.312 -1.107 -0.702 -0.506 -0.285 0.115 1.001 2200
## e.y[24] -0.976 0.273 -1.504 -1.157 -0.980 -0.792 -0.429 1.001 3000
## e.y[25] 1.324 0.293 0.744 1.129 1.320 1.521 1.906 1.001 3000
## e.y[26] -0.156 0.268 -0.678 -0.336 -0.164 0.027 0.387 1.001 3000
## e.y[27] -0.216 0.268 -0.737 -0.395 -0.224 -0.032 0.323 1.001 3000
## e.y[28] -0.903 0.401 -1.698 -1.165 -0.902 -0.633 -0.110 1.002 1400
## e.y[29] -2.492 0.302 -3.093 -2.692 -2.492 -2.298 -1.883 1.001 2900
## e.y[30] -0.781 0.422 -1.616 -1.058 -0.779 -0.501 0.060 1.002 1300
## e.y[31] 0.706 0.535 -0.331 0.343 0.699 1.076 1.730 1.002 1000
## e.y[32] -0.281 0.568 -1.392 -0.664 -0.285 0.096 0.854 1.003 940
## e.y[33] 0.167 0.596 -1.035 -0.217 0.157 0.556 1.346 1.003 950
## e.y[34] -0.609 0.460 -1.534 -0.913 -0.613 -0.311 0.309 1.002 1200
## e.y[35] -1.345 0.441 -2.224 -1.635 -1.346 -1.055 -0.474 1.002 1200
## e.y[36] 0.109 0.575 -1.059 -0.263 0.097 0.485 1.237 1.003 980
## e.y[37] 1.667 0.607 0.453 1.263 1.668 2.073 2.839 1.003 1000
## e.y[38] 1.427 0.611 0.225 1.022 1.430 1.831 2.609 1.002 1000
## e.y[39] 2.304 0.515 1.280 1.962 2.302 2.644 3.311 1.002 1200
## e.y[40] 1.958 0.435 1.114 1.666 1.956 2.248 2.820 1.002 1200
## e.y[41] 2.143 0.308 1.531 1.939 2.137 2.349 2.750 1.001 3000
## e.y[42] 3.197 0.352 2.506 2.966 3.199 3.433 3.903 1.001 3000
## e.y[43] 3.023 0.396 2.235 2.765 3.024 3.286 3.800 1.002 1400
## e.y[44] 0.294 0.476 -0.653 -0.029 0.306 0.600 1.231 1.001 2800
## e.y[45] -1.044 0.448 -1.935 -1.346 -1.031 -0.754 -0.168 1.001 3000
## e.y[46] -0.510 0.567 -1.647 -0.880 -0.506 -0.130 0.612 1.002 1100
## e.y[47] -0.666 0.558 -1.792 -1.026 -0.657 -0.298 0.452 1.003 1100
## e.y[48] -0.680 0.545 -1.759 -1.038 -0.680 -0.321 0.378 1.002 1100
## e.y[49] -0.671 0.322 -1.295 -0.883 -0.677 -0.451 -0.032 1.002 1900
## e.y[50] -0.698 0.294 -1.281 -0.894 -0.699 -0.498 -0.129 1.001 3000
## e.y[51] 2.185 0.442 1.281 1.898 2.184 2.484 3.063 1.001 3000
## e.y[52] -0.592 0.481 -1.535 -0.917 -0.576 -0.279 0.361 1.001 2500
## e.y[53] -3.643 0.604 -4.846 -4.052 -3.635 -3.233 -2.462 1.001 3000
## e.y[54] -0.398 0.459 -1.300 -0.707 -0.388 -0.096 0.518 1.001 2400
## sigma.y 1.924 0.196 1.589 1.785 1.908 2.040 2.363 1.001 3000
## deviance 222.200 2.935 218.548 220.085 221.490 223.655 229.525 1.001 3000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 4.3 and DIC = 226.5
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R² for pooling
SP_pooling <- jags_model_results$BUGSoutput$sims.list # Extract simulations from the results
r_y_pooling <- 1 - mean(apply(SP_pooling$e.y, 1, var)) / var(data$ICA) # Calculate R² for pooling
cat(paste("R² for pooling model:", round(r_y_pooling, 4)), "\n")
## R² for pooling model: 0.3402
# Load the devtools library
library(devtools)
## Loading required package: usethis
# Install the CalvinBayes package from GitHub
devtools::install_github("rpruim/CalvinBayes")
## Skipping install of 'CalvinBayes' from a github remote, the SHA1 (dbae67a3) has not changed since last install.
## Use `force = TRUE` to force installation
library(CalvinBayes)
## Loading required package: ggformula
## Loading required package: scales
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
## The following object is masked from 'package:arm':
##
## rescale
## Loading required package: ggridges
##
## New to ggformula? Try the tutorials:
## learnr::run_tutorial("introduction", package = "ggformula")
## learnr::run_tutorial("refining", package = "ggformula")
## Loading required package: bayesplot
## This is bayesplot version 1.11.1
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
##
## Attaching package: 'CalvinBayes'
## The following object is masked from 'package:bayesplot':
##
## rhat
## The following object is masked from 'package:datasets':
##
## HairEyeColor
# Define the parameters to analyze
parameters <- c("beta0", "beta1", "beta2")
# Generate and display diagnostic plots
for (param in parameters) {
diag_mcmc(as.mcmc(jags_model_results), parName = param)
}
Modelo 2: No pooling
# Define the JAGS model with factors 'epoca' and 'parte'
jags_model_code <- "
model {
for (i in 1:n) {
y[i] ~ dnorm(mu[i], tau.y) # Likelihood: y follows a normal distribution
mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i] + alpha_epoca[epoca[i]] + alpha_parte[parte[i]] # Linear predictor
e.y[i] <- y[i] - mu[i] # Residuals
}
beta0 ~ dnorm(0, 0.001) # Prior for intercept
beta1 ~ dnorm(0, 0.001) # Prior for slope of x1
beta2 ~ dnorm(0, 0.001) # Prior for slope of x2
# Priors for random effects associated with 'epoca'
for (j in 1:J_epoca) {
alpha_epoca[j] ~ dnorm(0, 0.001)
}
# Priors for random effects associated with 'parte'
for (k in 1:K_parte) {
alpha_parte[k] ~ dnorm(0, 0.001)
}
tau.y <- pow(sigma.y, -2) # Precision parameter
sigma.y ~ dunif(0, 100) # Prior for standard deviation
}"
writeLines(jags_model_code, con = "modelo_jags.bug") # Save model code to file
# Prepare data for JAGS
data_jags <- list(
y = data$ICA, # Response variable
x1 = as.numeric(data$`Temp agua (°C)`), # Predictor variable 1
x2 = as.numeric(data$`P sol (mgPO/l)`), # Predictor variable 2
epoca = data$Epoca, # Factor for random effect 'epoca'
parte = data$Parte, # Factor for random effect 'parte'
n = nrow(data), # Number of observations
J_epoca = max(as.numeric(as.factor(data$Epoca))), # Number of levels in 'epoca'
K_parte = max(as.numeric(as.factor(data$Parte))) # Number of levels in 'parte'
)
# Run the JAGS model
jags_model_results <- jags(
data = data_jags,
inits = function() list(beta0 = rnorm(1), beta1 = rnorm(1), sigma.y = runif(1, 0, 100)), # Initial values
parameters.to.save = c("beta0", "beta1", "beta2", "sigma.y", "alpha_epoca", "alpha_parte", "e.y"), # Parameters to monitor
model.file = "modelo_jags.bug", # Model file path
n.iter = 3000, # Total iterations
n.burnin = 1000, # Burn-in period
n.thin = 2, # Thinning interval
n.chains = 3 # Number of chains
)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 54
## Unobserved stochastic nodes: 16
## Total graph size: 470
##
## Initializing model
print(jags_model_results) # Print results
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
## 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 2
## n.sims = 3000 iterations saved. Running time = 0.202 secs
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat
## alpha_epoca[1] -2.457 15.533 -33.474 -12.606 -2.235 7.797 27.404 1.001
## alpha_epoca[2] -0.883 15.561 -31.928 -11.051 -0.824 9.363 28.799 1.001
## alpha_epoca[3] -3.840 15.538 -35.005 -13.997 -3.820 6.402 25.936 1.001
## alpha_parte[1] -0.767 10.048 -20.088 -7.515 -0.758 5.906 18.805 1.001
## alpha_parte[2] -1.934 10.051 -21.416 -8.696 -1.937 4.620 17.586 1.001
## alpha_parte[3] -1.228 10.034 -20.506 -7.941 -1.249 5.404 18.385 1.001
## alpha_parte[4] 0.112 10.028 -19.116 -6.568 0.070 6.781 19.895 1.001
## alpha_parte[5] -0.179 10.038 -19.511 -6.928 -0.184 6.533 19.297 1.001
## alpha_parte[6] -0.235 10.042 -19.327 -6.948 -0.220 6.405 19.426 1.001
## alpha_parte[7] -1.356 9.994 -20.506 -8.043 -1.259 5.194 18.132 1.001
## alpha_parte[8] -2.003 10.012 -21.384 -8.763 -1.946 4.549 17.670 1.001
## alpha_parte[9] -1.731 10.050 -20.996 -8.408 -1.744 4.843 17.836 1.001
## beta0 -6.387 17.684 -41.276 -18.250 -6.339 5.333 28.065 1.001
## beta1 0.593 0.183 0.235 0.471 0.589 0.718 0.951 1.001
## beta2 3.004 0.668 1.687 2.556 2.997 3.434 4.323 1.000
## e.y[1] -0.935 0.845 -2.608 -1.499 -0.933 -0.375 0.714 1.001
## e.y[2] -0.077 0.839 -1.669 -0.644 -0.076 0.462 1.599 1.002
## e.y[3] -1.115 0.843 -2.716 -1.681 -1.128 -0.556 0.553 1.001
## e.y[4] -0.883 0.836 -2.498 -1.465 -0.878 -0.336 0.768 1.001
## e.y[5] -1.651 0.821 -3.332 -2.194 -1.643 -1.097 -0.080 1.002
## e.y[6] -3.405 0.758 -4.906 -3.912 -3.398 -2.896 -1.943 1.001
## e.y[7] -0.555 0.799 -2.136 -1.094 -0.552 -0.038 1.042 1.001
## e.y[8] 3.205 0.721 1.816 2.730 3.201 3.685 4.594 1.007
## e.y[9] 3.203 0.754 1.707 2.699 3.213 3.692 4.712 1.003
## e.y[10] 0.981 0.841 -0.610 0.417 0.973 1.551 2.657 1.001
## e.y[11] 1.450 0.843 -0.249 0.895 1.452 2.003 3.086 1.002
## e.y[12] -0.672 0.816 -2.268 -1.200 -0.689 -0.138 0.945 1.001
## e.y[13] 0.679 0.729 -0.745 0.190 0.671 1.155 2.128 1.001
## e.y[14] -0.794 0.738 -2.278 -1.269 -0.790 -0.298 0.598 1.002
## e.y[15] -2.184 0.786 -3.748 -2.694 -2.174 -1.671 -0.644 1.002
## e.y[16] -1.927 0.796 -3.520 -2.462 -1.912 -1.400 -0.355 1.002
## e.y[17] -2.056 0.753 -3.556 -2.533 -2.064 -1.568 -0.568 1.001
## e.y[18] -1.939 0.787 -3.487 -2.450 -1.970 -1.386 -0.401 1.001
## e.y[19] -0.386 0.821 -1.985 -0.942 -0.389 0.171 1.214 1.001
## e.y[20] -0.392 0.837 -2.051 -0.937 -0.386 0.174 1.228 1.001
## e.y[21] -0.099 0.781 -1.600 -0.637 -0.099 0.425 1.450 1.001
## e.y[22] 0.023 0.732 -1.412 -0.478 0.021 0.489 1.493 1.001
## e.y[23] 0.314 0.749 -1.183 -0.167 0.318 0.806 1.789 1.001
## e.y[24] -0.432 0.771 -2.016 -0.927 -0.421 0.081 1.067 1.001
## e.y[25] 0.900 0.829 -0.723 0.348 0.893 1.457 2.503 1.002
## e.y[26] 0.593 0.785 -0.963 0.074 0.595 1.101 2.173 1.001
## e.y[27] -0.172 0.784 -1.696 -0.700 -0.185 0.368 1.411 1.001
## e.y[28] 0.239 0.833 -1.402 -0.324 0.242 0.778 1.923 1.001
## e.y[29] -1.469 0.818 -3.039 -1.995 -1.475 -0.926 0.181 1.001
## e.y[30] 0.782 0.819 -0.818 0.250 0.776 1.312 2.416 1.001
## e.y[31] 1.303 0.741 -0.134 0.808 1.278 1.792 2.771 1.001
## e.y[32] 0.723 0.803 -0.871 0.177 0.734 1.269 2.265 1.001
## e.y[33] 0.819 0.811 -0.881 0.302 0.830 1.351 2.393 1.001
## e.y[34] -0.437 0.809 -2.020 -0.971 -0.443 0.095 1.144 1.001
## e.y[35] 0.020 0.756 -1.471 -0.478 0.023 0.530 1.504 1.001
## e.y[36] 0.335 0.829 -1.267 -0.199 0.352 0.886 2.003 1.001
## e.y[37] 0.533 0.810 -1.057 -0.002 0.553 1.050 2.104 1.001
## e.y[38] 0.586 0.817 -1.085 0.052 0.593 1.123 2.176 1.001
## e.y[39] 1.229 0.779 -0.319 0.721 1.236 1.729 2.726 1.001
## e.y[40] -0.212 0.776 -1.684 -0.734 -0.213 0.275 1.354 1.001
## e.y[41] 1.273 0.743 -0.207 0.779 1.279 1.757 2.663 1.001
## e.y[42] 2.553 0.759 1.054 2.040 2.566 3.067 4.063 1.011
## e.y[43] 2.071 0.890 0.370 1.471 2.047 2.654 3.903 1.002
## e.y[44] 0.327 0.751 -1.146 -0.177 0.321 0.825 1.807 1.001
## e.y[45] -1.975 0.756 -3.500 -2.468 -1.971 -1.487 -0.477 1.001
## e.y[46] -0.390 0.747 -1.890 -0.865 -0.373 0.094 1.041 1.001
## e.y[47] -0.224 0.787 -1.797 -0.746 -0.222 0.302 1.319 1.001
## e.y[48] -0.091 0.801 -1.650 -0.628 -0.090 0.423 1.529 1.001
## e.y[49] -0.824 0.753 -2.269 -1.352 -0.827 -0.328 0.698 1.001
## e.y[50] 0.194 0.774 -1.279 -0.335 0.186 0.728 1.712 1.001
## e.y[51] 2.630 0.820 1.020 2.111 2.651 3.163 4.206 1.001
## e.y[52] -0.088 0.784 -1.625 -0.605 -0.101 0.430 1.476 1.001
## e.y[53] -2.140 0.800 -3.672 -2.673 -2.154 -1.626 -0.573 1.001
## e.y[54] 0.582 0.784 -0.963 0.068 0.567 1.110 2.123 1.001
## sigma.y 1.602 0.184 1.287 1.472 1.582 1.714 2.013 1.001
## deviance 202.299 6.277 192.411 197.760 201.505 205.842 216.747 1.001
## n.eff
## alpha_epoca[1] 3000
## alpha_epoca[2] 3000
## alpha_epoca[3] 3000
## alpha_parte[1] 3000
## alpha_parte[2] 3000
## alpha_parte[3] 3000
## alpha_parte[4] 3000
## alpha_parte[5] 3000
## alpha_parte[6] 3000
## alpha_parte[7] 3000
## alpha_parte[8] 3000
## alpha_parte[9] 3000
## beta0 3000
## beta1 3000
## beta2 3000
## e.y[1] 3000
## e.y[2] 1700
## e.y[3] 3000
## e.y[4] 3000
## e.y[5] 1900
## e.y[6] 3000
## e.y[7] 2800
## e.y[8] 3000
## e.y[9] 3000
## e.y[10] 3000
## e.y[11] 1700
## e.y[12] 2100
## e.y[13] 3000
## e.y[14] 1800
## e.y[15] 1400
## e.y[16] 1600
## e.y[17] 3000
## e.y[18] 3000
## e.y[19] 3000
## e.y[20] 3000
## e.y[21] 2200
## e.y[22] 3000
## e.y[23] 3000
## e.y[24] 2400
## e.y[25] 1600
## e.y[26] 3000
## e.y[27] 3000
## e.y[28] 3000
## e.y[29] 3000
## e.y[30] 3000
## e.y[31] 3000
## e.y[32] 3000
## e.y[33] 3000
## e.y[34] 2200
## e.y[35] 3000
## e.y[36] 3000
## e.y[37] 3000
## e.y[38] 2500
## e.y[39] 2000
## e.y[40] 3000
## e.y[41] 2700
## e.y[42] 1100
## e.y[43] 1000
## e.y[44] 3000
## e.y[45] 3000
## e.y[46] 3000
## e.y[47] 3000
## e.y[48] 2100
## e.y[49] 3000
## e.y[50] 3000
## e.y[51] 3000
## e.y[52] 2800
## e.y[53] 3000
## e.y[54] 3000
## sigma.y 2900
## deviance 2100
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 19.7 and DIC = 222.0
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R² for non-pooling model
SP_npooling <- jags_model_results$BUGSoutput$sims.list # Extract simulations from the results
r_y_npooling <- 1 - mean(apply(SP_npooling$e.y, 1, var)) / var(data$ICA) # Calculate R² for non-pooling model
cat(paste("R² for pooling model:", round(r_y_npooling, 4)), "\n")
## R² for pooling model: 0.5413
# Define the parameters to analyze
parameters <- c("beta0", "beta1", "beta2", "alpha_epoca[1]", "alpha_epoca[2]", "alpha_epoca[3]", "alpha_parte[1]", "alpha_parte[2]", "alpha_parte[3]", "alpha_parte[4]", "alpha_parte[5]", "alpha_parte[6]", "alpha_parte[7]", "alpha_parte[8]", "alpha_parte[9]")
# Generate and display diagnostic plots
for (param in parameters) {
diag_mcmc(as.mcmc(jags_model_results), parName = param)
}
Modelo 3: Pooling parcial multinivel no anidado e intercepto
variable
# Define the JAGS model
jags_model_code <- "
model {
for (i in 1:n) {
# Likelihood
y[i] ~ dnorm(mu[i], tau.y)
mu[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i] + alpha_epoca[epoca[i]] + alpha_parte[parte[i]]
e.y[i] <- y[i] - mu[i] # Residuals
}
# Priors for fixed effects
beta0 ~ dnorm(0, 0.001)
beta1 ~ dnorm(0, 0.001)
beta2 ~ dnorm(0, 0.001)
# Priors for random effects
for (j in 1:J_epoca) {
alpha_epoca[j] ~ dnorm(0, tau.epoca)
e.e[j] <- alpha_epoca[j] - beta0 # Residuals for epoch
}
for (k in 1:K_parte) {
alpha_parte[k] ~ dnorm(0, tau.parte)
e.p[k] <- alpha_parte[k] - beta0 # Residuals for part
}
# Priors for precisions
tau.y <- pow(sigma.y, -2)
tau.epoca <- pow(sigma.epoca, -2)
tau.parte <- pow(sigma.parte, -2)
sigma.y ~ dunif(0, 100) # Uniform prior for sigma.y
sigma.epoca ~ dunif(0, 100) # Uniform prior for sigma.epoca
sigma.parte ~ dunif(0, 100) # Uniform prior for sigma.parte
}
"
# Specify the path to save the JAGS model file
jags_model_path <- "modelo_jags.bug"
# Write the model code to a file
writeLines(jags_model_code, con = jags_model_path)
# Create a list of data for JAGS
data_jags <- list(
y = data$ICA, # Dependent variable
x1 = as.numeric(data$`Temp agua (°C)`), # Convert to numeric vector
x2 = as.numeric(data$`P sol (mgPO/l)`),
epoca = data$Epoca, # Epoch index
parte = data$Parte, # Part index
n = nrow(data), # Number of observations
J_epoca = max(as.numeric(as.factor(data$Epoca))), # Maximum epoch index
K_parte = max(as.numeric(as.factor(data$Parte))) # Number of levels in part
)
# Initial values function for JAGS model parameters
inits <- function() {
list(
beta0 = rnorm(1),
beta1 = rnorm(1),
alpha_epoca = rnorm(data_jags$J_epoca),
alpha_parte = rnorm(data_jags$K_parte),
sigma.y = runif(1, 0, 100),
sigma.epoca = runif(1, 0, 100),
sigma.parte = runif(1, 0, 100)
)
}
# Parameters to monitor during the JAGS model run
params <- c("beta0", "beta1", "beta2", "alpha_epoca", "alpha_parte", "sigma.y", "sigma.epoca", "sigma.parte", "e.e", "e.p", "e.y")
# Load and run the JAGS model
jags_model_results <- jags(
data = data_jags,
inits = inits,
parameters.to.save = params,
model.file = jags_model_path,
n.iter = 3000, # Total number of iterations
n.burnin = 1000, # Number of burn-in iterations
n.thin = 2, # Thinning interval (sample every second iteration)
n.chains = 3 # Number of chains to run
)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 54
## Unobserved stochastic nodes: 18
## Total graph size: 486
##
## Initializing model
# Print summary of the results from the JAGS model
print(jags_model_results)
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
## 3 chains, each with 3000 iterations (first 1000 discarded), n.thin = 2
## n.sims = 3000 iterations saved. Running time = 0.216 secs
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat
## alpha_epoca[1] -0.245 5.461 -11.679 -1.335 -0.133 0.996 9.914 1.008
## alpha_epoca[2] 1.368 5.457 -9.945 0.274 1.425 2.609 11.248 1.008
## alpha_epoca[3] -1.433 5.465 -12.882 -2.529 -1.288 -0.210 8.564 1.008
## alpha_parte[1] 0.091 0.341 -0.570 -0.072 0.033 0.243 0.898 1.004
## alpha_parte[2] -0.156 0.378 -1.129 -0.305 -0.066 0.035 0.463 1.001
## alpha_parte[3] 0.014 0.343 -0.715 -0.121 0.004 0.162 0.760 1.004
## alpha_parte[4] 0.165 0.383 -0.426 -0.034 0.069 0.317 1.132 1.004
## alpha_parte[5] 0.105 0.358 -0.519 -0.069 0.032 0.245 0.985 1.001
## alpha_parte[6] 0.106 0.356 -0.531 -0.060 0.035 0.239 1.031 1.002
## alpha_parte[7] -0.045 0.338 -0.845 -0.185 -0.010 0.106 0.663 1.001
## alpha_parte[8] -0.205 0.392 -1.256 -0.368 -0.088 0.015 0.371 1.003
## alpha_parte[9] -0.113 0.350 -0.983 -0.254 -0.035 0.059 0.496 1.002
## beta0 -5.879 6.142 -17.104 -8.454 -5.957 -3.318 6.037 1.007
## beta1 0.447 0.132 0.193 0.359 0.444 0.531 0.710 1.001
## beta2 2.304 0.577 1.211 1.908 2.291 2.670 3.480 1.001
## e.e[1] 5.634 11.261 -17.174 2.414 5.835 9.102 25.649 1.008
## e.e[2] 7.247 11.235 -15.456 4.017 7.439 10.683 27.616 1.008
## e.e[3] 4.446 11.278 -18.915 1.202 4.736 7.959 24.228 1.008
## e.p[1] 5.970 6.138 -5.847 3.371 6.055 8.560 17.147 1.007
## e.p[2] 5.723 6.107 -6.344 3.240 5.758 8.286 16.745 1.007
## e.p[3] 5.893 6.127 -5.989 3.323 6.003 8.513 16.825 1.007
## e.p[4] 6.044 6.202 -5.996 3.421 6.084 8.639 17.437 1.007
## e.p[5] 5.984 6.203 -6.049 3.287 6.024 8.584 17.369 1.007
## e.p[6] 5.985 6.206 -6.084 3.314 6.068 8.598 17.307 1.007
## e.p[7] 5.834 6.166 -6.173 3.266 5.884 8.440 16.995 1.007
## e.p[8] 5.674 6.150 -6.435 3.111 5.743 8.321 16.865 1.007
## e.p[9] 5.766 6.137 -6.208 3.220 5.778 8.354 16.800 1.007
## e.y[1] -0.969 0.672 -2.337 -1.410 -0.960 -0.521 0.319 1.001
## e.y[2] -0.236 0.625 -1.462 -0.656 -0.231 0.176 0.988 1.001
## e.y[3] -1.078 0.634 -2.329 -1.474 -1.072 -0.659 0.097 1.001
## e.y[4] -1.598 0.554 -2.687 -1.955 -1.600 -1.255 -0.497 1.001
## e.y[5] -2.751 0.563 -3.770 -3.133 -2.785 -2.401 -1.527 1.002
## e.y[6] -3.961 0.490 -4.940 -4.275 -3.975 -3.653 -2.949 1.002
## e.y[7] -0.060 0.600 -1.302 -0.438 -0.040 0.329 1.072 1.001
## e.y[8] 2.634 0.528 1.599 2.289 2.621 2.971 3.736 1.002
## e.y[9] 3.241 0.544 2.127 2.891 3.249 3.597 4.301 1.001
## e.y[10] 1.854 0.596 0.602 1.459 1.879 2.264 2.966 1.001
## e.y[11] 2.049 0.593 0.858 1.667 2.059 2.449 3.202 1.001
## e.y[12] -0.086 0.583 -1.279 -0.467 -0.069 0.301 1.046 1.001
## e.y[13] 0.483 0.504 -0.481 0.146 0.486 0.805 1.476 1.001
## e.y[14] -1.446 0.527 -2.435 -1.802 -1.457 -1.113 -0.329 1.003
## e.y[15] -2.628 0.512 -3.591 -2.977 -2.631 -2.309 -1.583 1.001
## e.y[16] -1.174 0.589 -2.371 -1.558 -1.141 -0.786 -0.044 1.001
## e.y[17] -2.151 0.595 -3.275 -2.548 -2.156 -1.780 -0.961 1.001
## e.y[18] -1.454 0.606 -2.662 -1.849 -1.450 -1.061 -0.262 1.001
## e.y[19] 0.512 0.564 -0.680 0.165 0.528 0.894 1.583 1.002
## e.y[20] 0.349 0.538 -0.778 0.009 0.367 0.720 1.361 1.001
## e.y[21] 0.527 0.542 -0.617 0.175 0.546 0.890 1.547 1.002
## e.y[22] -0.322 0.518 -1.335 -0.666 -0.321 0.021 0.703 1.001
## e.y[23] -0.431 0.528 -1.390 -0.793 -0.447 -0.094 0.697 1.001
## e.y[24] -0.881 0.496 -1.831 -1.217 -0.883 -0.570 0.123 1.001
## e.y[25] 1.139 0.511 0.062 0.816 1.147 1.493 2.116 1.001
## e.y[26] 0.028 0.504 -0.949 -0.313 0.016 0.347 1.058 1.001
## e.y[27] -0.187 0.492 -1.176 -0.502 -0.179 0.139 0.743 1.002
## e.y[28] 0.662 0.600 -0.535 0.267 0.680 1.068 1.817 1.001
## e.y[29] -1.042 0.543 -2.141 -1.380 -1.026 -0.685 -0.002 1.001
## e.y[30] 0.871 0.613 -0.326 0.465 0.866 1.279 2.071 1.001
## e.y[31] 1.257 0.531 0.194 0.911 1.267 1.614 2.282 1.001
## e.y[32] 0.486 0.572 -0.645 0.102 0.490 0.866 1.606 1.002
## e.y[33] 0.915 0.596 -0.290 0.510 0.924 1.309 2.076 1.001
## e.y[34] 0.256 0.579 -0.931 -0.102 0.275 0.653 1.332 1.001
## e.y[35] -0.277 0.528 -1.286 -0.630 -0.280 0.060 0.798 1.002
## e.y[36] 0.803 0.615 -0.450 0.413 0.817 1.222 1.964 1.001
## e.y[37] 0.819 0.620 -0.397 0.396 0.828 1.238 2.011 1.001
## e.y[38] 0.701 0.628 -0.494 0.273 0.697 1.118 1.944 1.001
## e.y[39] 1.381 0.573 0.249 0.995 1.381 1.755 2.498 1.001
## e.y[40] 0.039 0.597 -1.162 -0.348 0.043 0.427 1.184 1.001
## e.y[41] 0.847 0.517 -0.119 0.501 0.835 1.178 1.922 1.002
## e.y[42] 2.182 0.500 1.241 1.848 2.167 2.500 3.198 1.002
## e.y[43] 1.867 0.553 0.812 1.508 1.854 2.227 2.978 1.001
## e.y[44] -0.333 0.573 -1.373 -0.720 -0.354 0.023 0.862 1.002
## e.y[45] -1.959 0.518 -2.974 -2.310 -1.956 -1.623 -0.949 1.002
## e.y[46] -0.195 0.579 -1.322 -0.587 -0.184 0.183 0.953 1.001
## e.y[47] -0.228 0.587 -1.410 -0.627 -0.226 0.168 0.948 1.001
## e.y[48] 0.077 0.592 -1.096 -0.321 0.088 0.468 1.257 1.001
## e.y[49] -0.723 0.510 -1.726 -1.050 -0.720 -0.378 0.245 1.001
## e.y[50] -0.274 0.516 -1.228 -0.619 -0.300 0.048 0.834 1.001
## e.y[51] 2.655 0.571 1.523 2.284 2.643 3.021 3.781 1.001
## e.y[52] 0.169 0.564 -0.961 -0.207 0.180 0.530 1.308 1.001
## e.y[53] -2.508 0.663 -3.766 -2.945 -2.519 -2.084 -1.208 1.001
## e.y[54] 0.403 0.553 -0.683 0.038 0.407 0.758 1.498 1.001
## sigma.epoca 6.501 10.346 0.779 1.694 2.938 6.174 38.832 1.001
## sigma.parte 0.350 0.287 0.004 0.134 0.288 0.497 1.070 1.161
## sigma.y 1.570 0.166 1.292 1.451 1.556 1.676 1.927 1.001
## deviance 200.126 4.394 192.984 197.144 199.519 202.534 210.752 1.002
## n.eff
## alpha_epoca[1] 3000
## alpha_epoca[2] 3000
## alpha_epoca[3] 3000
## alpha_parte[1] 1800
## alpha_parte[2] 3000
## alpha_parte[3] 1100
## alpha_parte[4] 780
## alpha_parte[5] 3000
## alpha_parte[6] 1100
## alpha_parte[7] 3000
## alpha_parte[8] 1500
## alpha_parte[9] 3000
## beta0 3000
## beta1 3000
## beta2 3000
## e.e[1] 3000
## e.e[2] 3000
## e.e[3] 3000
## e.p[1] 2700
## e.p[2] 3000
## e.p[3] 2700
## e.p[4] 2800
## e.p[5] 2900
## e.p[6] 3000
## e.p[7] 3000
## e.p[8] 3000
## e.p[9] 2700
## e.y[1] 3000
## e.y[2] 3000
## e.y[3] 3000
## e.y[4] 3000
## e.y[5] 1000
## e.y[6] 1800
## e.y[7] 2900
## e.y[8] 1200
## e.y[9] 2800
## e.y[10] 3000
## e.y[11] 3000
## e.y[12] 3000
## e.y[13] 3000
## e.y[14] 930
## e.y[15] 3000
## e.y[16] 3000
## e.y[17] 3000
## e.y[18] 3000
## e.y[19] 1800
## e.y[20] 3000
## e.y[21] 2200
## e.y[22] 3000
## e.y[23] 2200
## e.y[24] 3000
## e.y[25] 3000
## e.y[26] 3000
## e.y[27] 2300
## e.y[28] 3000
## e.y[29] 3000
## e.y[30] 3000
## e.y[31] 3000
## e.y[32] 1700
## e.y[33] 3000
## e.y[34] 3000
## e.y[35] 1900
## e.y[36] 3000
## e.y[37] 3000
## e.y[38] 3000
## e.y[39] 3000
## e.y[40] 3000
## e.y[41] 1100
## e.y[42] 3000
## e.y[43] 3000
## e.y[44] 1700
## e.y[45] 3000
## e.y[46] 3000
## e.y[47] 3000
## e.y[48] 3000
## e.y[49] 3000
## e.y[50] 2900
## e.y[51] 3000
## e.y[52] 3000
## e.y[53] 3000
## e.y[54] 3000
## sigma.epoca 3000
## sigma.parte 37
## sigma.y 3000
## deviance 1100
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 9.6 and DIC = 209.8
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R²
SP_model3 <- jags_model_results$BUGSoutput$sims.list # Extract simulations from the results
r_y_model3 <- 1 - mean(apply(SP_model3$e.y, 1, var)) / var(data$ICA) # Calculate R²
cat(paste("R² for model 3:", round(r_y_model3, 4)), "\n")
## R² for model 3: 0.561
# Extract posterior samples and calculate omega
posterior_samples <- jags_model_results$BUGSoutput$sims.list
omega_model3 <- pmin((apply(posterior_samples$e.y, 2, sd) / mean(posterior_samples$sigma.y))^2, 1)
cat("Omega:\n")
## Omega:
plot(
omega_model3,
type = "l", col = "blue", lwd = 2,
main = "", xlab = "Index", ylab = "Omega"
)
# Define the parameters to analyze
parameters <- c("beta0", "beta1", "beta2", "alpha_epoca[1]", "alpha_epoca[2]", "alpha_epoca[3]", "alpha_parte[1]", "alpha_parte[2]", "alpha_parte[3]", "alpha_parte[4]", "alpha_parte[5]", "alpha_parte[6]", "alpha_parte[7]", "alpha_parte[8]", "alpha_parte[9]")
# Generate and display diagnostic plots
for (param in parameters) {
diag_mcmc(as.mcmc(jags_model_results), parName = param)
}
Modelo 4: Pooling parcial multinivel con efectos aleatorios a nivel
de pendientes
# Define the JAGS model
jags_model_code <- "
model {
for (i in 1:n) {
# Likelihood
y[i] ~ dnorm(mu[i], tau.y)
mu[i] <- beta0 + beta1[epoca[i]] * x1[i] + beta2[epoca[i]] * x2[i]
e.y[i] <- y[i] - mu[i] # Residuals
}
# Priors for fixed effects
beta0 ~ dnorm(0, 0.001)
# Priors for precision of the likelihood
tau.y <- pow(sigma.y, -2)
sigma.y ~ dunif(0, 100)
# Priors for random effects
for (j in 1:J_epoca) {
beta1[j] ~ dnorm(m1, t1)
beta2[j] ~ dnorm(m2, t2)
}
# Hyperparameters for random effects
m1 ~ dnorm(0, 0.0001)
m2 ~ dnorm(0, 0.0001)
t1 <- pow(s1, -2)
s1 ~ dunif(0, 100)
t2 <- pow(s2, -2)
s2 ~ dunif(0, 100)
}
"
# Specify the path to save the JAGS model file
jags_model_path <- "modelo_jags.bug"
# Write the model code to a file
writeLines(jags_model_code, con = jags_model_path)
# Create a list of data for JAGS
data_jags <- list(
y = data$ICA, # Dependent variable
x1 = as.numeric(data$`Temp agua (°C)`), # Independent variable 1
x2 = as.numeric(data$`P sol (mgPO/l)`), # Independent variable 2
epoca = data$Epoca, # Epoch index
n = nrow(data), # Number of observations
J_epoca = max(as.numeric(as.factor(data$Epoca))) # Maximum epoch index
)
# Initial values function for JAGS model parameters
inits <- function() {
list(
sigma.y = runif(1, 0, 100), # Initial value for sigma.y
s1 = runif(1, 0, 100), # Initial value for s1
s2 = runif(1, 0, 100) # Initial value for s2
)
}
# Parameters to monitor during the JAGS model run
params <- c("beta0", "beta1", "beta2", "sigma.y", "s1", "s2", "m1", "m2", "e.y")
# Load and run the JAGS model
jags_model_results <- jags(
data = data_jags,
inits = inits,
parameters.to.save = params,
model.file = jags_model_path,
n.iter = 10000, # Total number of iterations
n.burnin = 3000, # Number of burn-in iterations
n.thin = 4, # Thinning interval (sample every second iteration)
n.chains = 3 # Number of chains to run
)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 54
## Unobserved stochastic nodes: 12
## Total graph size: 424
##
## Initializing model
# Print summary of the results from the JAGS model
print(jags_model_results)
## Inference for Bugs model at "modelo_jags.bug", fit using jags,
## 3 chains, each with 10000 iterations (first 3000 discarded), n.thin = 4
## n.sims = 5250 iterations saved. Running time = 0.342 secs
## mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
## beta0 -6.031 2.382 -10.305 -7.858 -6.093 -4.344 -1.317 1.027 79
## beta1[1] 0.439 0.114 0.212 0.357 0.442 0.523 0.652 1.025 93
## beta1[2] 0.526 0.112 0.307 0.445 0.529 0.610 0.730 1.026 84
## beta1[3] 0.399 0.124 0.158 0.308 0.404 0.492 0.624 1.023 90
## beta2[1] 2.405 0.856 0.843 1.841 2.352 2.926 4.245 1.001 5200
## beta2[2] 1.895 0.809 0.265 1.375 1.923 2.417 3.468 1.001 2800
## beta2[3] 2.156 0.771 0.656 1.654 2.147 2.654 3.727 1.001 3300
## e.y[1] -1.164 0.591 -2.335 -1.570 -1.160 -0.756 -0.029 1.017 130
## e.y[2] -0.459 0.549 -1.549 -0.830 -0.451 -0.085 0.598 1.014 150
## e.y[3] -1.250 0.556 -2.326 -1.631 -1.252 -0.864 -0.188 1.019 110
## e.y[4] -1.849 0.440 -2.702 -2.141 -1.852 -1.555 -0.990 1.012 180
## e.y[5] -3.128 0.405 -3.918 -3.399 -3.128 -2.860 -2.324 1.007 300
## e.y[6] -4.125 0.425 -4.959 -4.412 -4.128 -3.840 -3.271 1.002 1600
## e.y[7] 0.077 0.577 -1.042 -0.327 0.077 0.466 1.229 1.009 250
## e.y[8] 2.476 0.492 1.517 2.141 2.475 2.806 3.451 1.005 430
## e.y[9] 3.277 0.533 2.235 2.909 3.278 3.637 4.337 1.007 310
## e.y[10] 2.380 0.504 1.403 2.041 2.373 2.716 3.373 1.010 230
## e.y[11] 2.538 0.523 1.533 2.187 2.532 2.892 3.569 1.011 200
## e.y[12] 0.380 0.504 -0.597 0.041 0.373 0.716 1.373 1.010 210
## e.y[13] 0.696 0.443 -0.155 0.398 0.691 0.993 1.569 1.004 570
## e.y[14] -1.409 0.438 -2.250 -1.701 -1.414 -1.116 -0.553 1.003 740
## e.y[15] -2.514 0.433 -3.342 -2.804 -2.519 -2.224 -1.669 1.003 1000
## e.y[16] -1.093 0.495 -2.044 -1.426 -1.099 -0.761 -0.106 1.005 430
## e.y[17] -2.356 0.526 -3.346 -2.709 -2.368 -1.997 -1.316 1.007 310
## e.y[18] -1.514 0.546 -2.547 -1.885 -1.525 -1.144 -0.424 1.008 260
## e.y[19] 0.634 0.470 -0.311 0.330 0.642 0.949 1.563 1.008 270
## e.y[20] 0.414 0.451 -0.480 0.123 0.418 0.716 1.291 1.006 400
## e.y[21] 0.590 0.466 -0.351 0.289 0.598 0.904 1.508 1.008 290
## e.y[22] -0.410 0.466 -1.351 -0.711 -0.402 -0.096 0.508 1.008 290
## e.y[23] -0.674 0.444 -1.562 -0.965 -0.670 -0.374 0.189 1.005 480
## e.y[24] -1.025 0.432 -1.889 -1.307 -1.023 -0.731 -0.199 1.002 1500
## e.y[25] 1.194 0.437 0.318 0.905 1.199 1.488 2.041 1.004 690
## e.y[26] -0.157 0.432 -1.016 -0.440 -0.157 0.138 0.664 1.001 2900
## e.y[27] -0.201 0.432 -1.067 -0.487 -0.202 0.096 0.618 1.001 3700
## e.y[28] 0.637 0.679 -0.654 0.188 0.636 1.088 1.987 1.001 5200
## e.y[29] -1.070 0.551 -2.135 -1.447 -1.076 -0.696 0.026 1.001 4700
## e.y[30] 0.774 0.703 -0.572 0.311 0.768 1.239 2.173 1.001 5200
## e.y[31] 1.179 0.568 0.024 0.807 1.185 1.545 2.286 1.003 1100
## e.y[32] 0.334 0.620 -0.906 -0.059 0.337 0.733 1.538 1.001 5200
## e.y[33] 0.901 0.644 -0.374 0.483 0.904 1.318 2.154 1.002 2400
## e.y[34] 0.417 0.515 -0.589 0.074 0.420 0.749 1.435 1.003 780
## e.y[35] -0.385 0.497 -1.365 -0.712 -0.382 -0.064 0.607 1.002 1300
## e.y[36] 0.919 0.621 -0.307 0.516 0.917 1.323 2.128 1.002 1600
## e.y[37] 0.823 0.748 -0.667 0.328 0.835 1.316 2.268 1.001 5200
## e.y[38] 0.612 0.757 -0.903 0.119 0.625 1.111 2.065 1.001 5200
## e.y[39] 1.382 0.634 0.129 0.962 1.389 1.797 2.617 1.001 5200
## e.y[40] 0.382 0.688 -0.977 -0.058 0.360 0.830 1.749 1.002 1900
## e.y[41] 0.820 0.453 -0.058 0.517 0.809 1.126 1.698 1.001 5200
## e.y[42] 2.076 0.411 1.275 1.798 2.072 2.342 2.889 1.002 2000
## e.y[43] 1.962 0.489 1.004 1.630 1.965 2.284 2.944 1.002 2300
## e.y[44] -0.640 0.518 -1.636 -0.992 -0.641 -0.301 0.394 1.004 640
## e.y[45] -2.046 0.472 -2.943 -2.361 -2.054 -1.740 -1.106 1.004 530
## e.y[46] 0.011 0.651 -1.231 -0.420 0.000 0.425 1.346 1.006 420
## e.y[47] -0.078 0.661 -1.334 -0.516 -0.096 0.340 1.281 1.004 550
## e.y[48] 0.247 0.719 -1.093 -0.244 0.229 0.714 1.744 1.001 5200
## e.y[49] -0.815 0.530 -1.871 -1.158 -0.807 -0.452 0.198 1.001 5200
## e.y[50] -0.496 0.440 -1.359 -0.788 -0.497 -0.193 0.356 1.001 2900
## e.y[51] 2.526 0.551 1.438 2.159 2.527 2.893 3.606 1.005 480
## e.y[52] 0.303 0.589 -0.840 -0.100 0.297 0.710 1.452 1.003 760
## e.y[53] -2.621 0.662 -3.916 -3.074 -2.630 -2.167 -1.315 1.008 290
## e.y[54] 0.459 0.573 -0.637 0.066 0.452 0.851 1.579 1.003 1000
## m1 0.460 0.838 -0.336 0.340 0.461 0.575 1.100 1.081 5200
## m2 2.207 4.589 -2.868 1.521 2.166 2.795 7.467 1.066 2700
## s1 0.389 1.156 0.026 0.078 0.140 0.303 2.484 1.004 770
## s2 3.023 7.570 0.047 0.432 1.043 2.516 19.919 1.001 4200
## sigma.y 1.630 0.171 1.340 1.506 1.615 1.734 2.008 1.001 5200
## deviance 204.454 4.014 198.692 201.449 203.766 206.802 213.816 1.001 5200
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule: pV = var(deviance)/2)
## pV = 8.1 and DIC = 212.5
## DIC is an estimate of expected predictive error (lower deviance is better).
# Calculate R²
SP_model4 <- jags_model_results$BUGSoutput$sims.list # Extract simulations from the results
r_y_model4 <- 1 - mean(apply(SP_model4$e.y , 1, var)) / var(data$ICA) # Calculate R²
cat(paste("R² for model 4:", round(r_y_model4, 4)), "\n")
## R² for model 4: 0.5247
# Extract posterior samples and calculate omega
posterior_samples <- jags_model_results$BUGSoutput$sims.list
omega_model4 <- pmin((apply(posterior_samples$e.y, 2, sd) / mean(posterior_samples$sigma.y))^2, 1)
cat("Omega:\n")
## Omega:
plot(
omega_model4,
type = "l", col = "blue", lwd = 2,
main = "", xlab = "Index", ylab = "Omega"
)
# Define the parameters to analyze
parameters <- c("beta0", "beta1[1]", "beta1[2]", "beta1[3]", "beta2[1]", "beta2[2]", "beta2[3]")
# Generate and display diagnostic plots
for (param in parameters) {
diag_mcmc(as.mcmc(jags_model_results), parName = param)
}
